Lattice study of a Janus interface 
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A lattice gas simulation of water between a hydrophobic plate and a hydrophilic plate (a Janus 
interface) shows large fluctuations in the number of liquid cells in contact with the hydrophobic 
plate, and a power spectrum similar to the experimental results that Zhang, Zhu, and Granick 
found [X. Y. Zhang, Y. X. Zhu, and S. Granick, Science 295, 663 (2002)] when measuring viscous 
response in a Janus system. Study of the spatial Fourier modes of the liquid- vapor interface suggests 
that interface fluctuations with length scales between approximately 1.5 and 20 nm cause the effects 
observed in the simulation. 



PACS numbers: 68.15.+e, 68.08.Bc 



I. INTRODUCTION 

The interactions between water and hydrophobic and 
hydrophilic surfaces are important in many areas of re- 
search, including such diverse fields as lubrication and 
protein folding. In order to study hydrophobic effects 
with as few additional complications as possible, the be- 
havior of water confined between two hydrophobic plates 
has been the subject of many experimental and theo- 
retical studies 0, IE IE S IE IE 0- More recently, experi- 
mentalists and theorists have also examined the behavior 
of water confined between one hydrophobic and one hy- 
drophilic plate[EEi, a so-called Janus interface, to learn 
how water reacts to the competing effects of the plates. 

While performing experiments to study the response 
of water in a Janus system to shear deformations, Zhang 
et.alp| found that the fluctuations about the mean value 
of the viscous response were typically 25-50% of the mean 
value, but ranged up to 100%. These fluctuations are ex- 
traordinarily large compared to the fluctuations observed 
when both plates are hydrophilic, and it is surprising that 
they do not average out over the plates, which are ap- 
proximately 10 (im on a side. When the time series of 
the viscous response was Fourier analyzed, it was found 
that the resulting power spectrum is independent of fre- 
quency at low frequencies, but at approximately / = .001 
Hz, the power spectrum crosses over to power-law depen- 
dence on frequency, with an apparent exponent near /~ 2 . 
For frequencies greater than / = .01 Hz, the power spec- 
trum again levels off and becomes seemingly independent 
of frequency. 

It has been predicted 01 that when water encounters a 
large hydrophobic solute, such as the hydrophobic plate 
in a Janus system, the water forms a semi-free liquid- 
vapor interface near the hydrophobic surface. The sur- 
face is unable to form hydrogen bonds with the water, 
and the formation of an interface allows the water to 
form a more energetically favorable hydrogen bonding 
network. Zhang et.al. suggested that their results are 
caused by fluctuations of such an interface. They further 
identified the time scale of ~ 10 3 /27r seconds, the inverse 
of the frequency at which power-law behavior begins in 
their power spectrum, as the lifetime of vapor or a vapor 



bubble. 

This paper presents the results of a simulation designed 
to test the hypothesis that fluctuations in the liquid- 
vapor interface are the source of the large fluctuations 
and power spectrum seen experimentally. We chose to 
study a lattice gas because it is a simple model system 
that can be used to examine the large length scale be- 
havior of a liquid between competing solvophobic and 
solvophilic plates, and can also describe interface fluctu- 
ations. After describing the details of the simulation in 
Section [HJ results are presented in Section ITTT1 and the 
paper is concluded with a brief discussion. 



II. 



SIMULATION DETAILS 



The system in the simulation is an L X L X R cubic 
lattice gas, where L is the side length of the hydrophobic 
and hydrophilic plates in the x and y directions in units 
of the lattice spacing a, and R is the distance between the 
plates, in the z direction, as shown in Figure^ L and R 
were chosen to be 128 and 8, respectively. Each cell i has 
an associated variable n, which labels the cell as liquid 
(rii = 1) or vapor (n, = 0). The cells around the edges 
of the system are constrained to have rii — 1 in order to 
ensure that there is no evaporation between the plates. 
The cells in the z = and z = 7 layers represent the 
hydrophobic and hydrophilic plates, respectively, and are 
constrained to have rij = 1 at all times so that liquid in 
the bulk always interacts with the plates. The remaining 
unconstrained cells in the bulk obey the Hamiltonian 



(3H = -Bey^njUj -f3aA^2m -/fos^n, 

<ij> i i n.n. A i n.n. B 

(i) 

where (3 = (ksT) , ks is Boltzmann's constant, < ij > 
indicates a sum over nearest neighbors, and i n.n. A and 
i n.n. B indicate sums over cells adjacent to plate A (the 
hydrophobic plate) or B (the hydrophilic plate). 

The parameters (3e and /3/z have the same values as 
in the simulation of Luzar and Leung 7], (3e — 1.26 and 
f3fi = — 6e/2 + 1.84 x 10~ 4 , chosen to match the surface 
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tension and isothermal compressibility of water at 300 
K. With these parameters, the lattice spacing a corre- 
sponds to .193 nm, so the plates are approximately 1.5 
nm apart, within the range studied by Zhang et.al. How- 
ever, the plates in the simulation have a side length of 
approximately 25 nm, which is significantly smaller than 
the experimental plates. The interactions with the plates 
are (3ub — /3e and (3a a — -l/3e, which gives a contact an- 
gle of 143° at the hydrophobic plate, comparable to the 
contact angle of approximately 120° in the experimental 
system. 

Zhang et.al. were careful to show that although they 
applied shear forces to the hydrophobic plate and mea- 
sured the response of the water in their experiments, the 
system was still in the linear response regime, so fluctua- 
tions present in the nonequilibrium experiment should be 
present at equilibrium as well. To make the simulation as 
simple as possible, we chose to study an equilibrium lat- 
tice gas propagated by Metropolis Monte Carlo, accord- 
ing to the Hamiltonian in equation Q for approximately 
140,000 passes through the lattice. The contact density, 
defined to be the number of liquid cells in contact with 
the hydrophobic plate, 



c(t)= 



(2) 



i n.n. A 



was calculated as a function of time, where time is mea- 
sured in units of Monte Carlo passes through the lattice. 
c(t) should be related to the viscosity of water near the 
hydrophobic plate, because the viscosity measured at the 
plate should increase with the amount of liquid (rather 
than vapor) in contact with the plate. 

The distance of the liquid-vapor interface h from the 
hydrophobic plate, as a function of (x, y) and time, was 
also calculated in the simulation. The interface in the 
(x, y) column is defined to be located at the average of 
the distances of the liquid cell closest to, and the vapor 
cell farthest from, the hydrophobic plate in that column, 
reminiscent of Weeks' derivation of the capillary wave 
model[ll[l2: 



Kx,y) = ^ 



max [Zi(l 



mm [ZiTii 
(x,y) 



(3) 



The value of h when liquid is in contact with the hy- 
drophobic plate is h = (0+l)/2 = 0.5, and increases in in- 
crements of 1 to its maximal value of h = (6 + 7)/2 = 6.5, 
for vapor in contact with the hydrophilic plate. 

The fact that the edges of the system are constrained 
to be liquid for all z leads to a boundary region of approx- 
imately 12 cells from each edge where the fluctuations of 
h are smaller than those in the middle of the lattice. Re- 
sults reported below for the 104 x 104 x 8 lattice resulting 
from ignoring the outer perimeter of 12 cells are virtu- 
ally indistinguishable from results with a boundary of 32 
cells, so the lattice seems to be large enough to study a 
Janus interface. 



Inside the boundary region, the time and space- 
averaged interface position is (h) = 2.0, which means 
that on average, the first one or two layers in contact 
with the hydrophobic plate are vapor. Because the plate 
is so significantly dewetted, changes in the amount of 
liquid in contact with the plate in our simulation are 
caused largely by fluctuations of the liquid-vapor inter- 
face, rather than by spontaneously-forming vapor bub- 
bles. To the extent that this lattice gas is an accu- 
rate description of the experimental Janus system, then, 
spontaneously-forming vapor bubbles do not seem to be 
important contributors to the results. 



III. 



RESULTS 



The contact density c(t) and interface h(x, y; t) were 
saved every 10 passes through the lattice during the simu- 
lation. A section of the time series for the contact density 
measured over the central 104 x 104 x 8 lattice is shown 
in Figure |2 The fluctuations in c(t) range up to approx- 
imately 25% of the mean, which is within the range seen 
experimentally. The power spectrum resulting from the 
time Fourier transform of c(t), S(f) = |c(/)| 2 , where 



c(/) = (At) ]T c{t)exp[2irft/t max ] 



(4) 



and At is the time between samples (10 passes in our 
simulation) is shown on a log-log plot in Figure The 
power is independent of frequency until approximately 
/ = .001 passes -1 , where it begins to drop off. The 
best- fit line to the data between / = .001 passes -1 and 
/ = .015 passes -1 has a slope of — 1.48±.04, less than the 
apparent power-law of / -2 seen experimentally. There 
seems to be a slight leveling off of the power as a func- 
tion of frequency at high /, but the data are so noisy 
that it is difficult to quantify. The contact density c(t) 
thus exhibits behavior surprisingly similar to the viscous 
response in the experiments of Zhang et.al. 

In order to try to understand the cause of the fluctu- 
ations and power spectrum of c(t), the dynamics of the 
liquid-vapor interface were studied. The characteristic 
length scales of fluctuations of the liquid- vapor interface 
relevant to fluctuations in the contact density can be in- 
ferred from an analysis of the spatial Fourier modes of the 
interface. The Fourier modes h(k x , k y ;t) of the interface 
can be defined with the equation 



JV 

h(k x , k y ;t) = (Ax) 2 M*> V, t)e- 2 ^ k * x+k ^ (5) 

x iV— o 

where Aa; is the spatial resolution (a in the simulation) 
and N is the number of cells in each dimension (104 in 
the simulation). 
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Figure 01 shows the time autocorrelation function of the 
k = (0, 0) spatial Fourier mode of the interface, 

F(t) = (Sh(0, 0; 0) 5h(0, 0; t)\ (6) 

where 5h(k x , k y ;t) = h(k x , k y ;t) — (h(k x ,k y ;t)) , as a 
function of time. F(t) decays on a time scale of approxi- 
mately 1000 passes, which is the inverse of the frequency 
near which power-law behavior begins in S(f). In the 
sense that fluctuations of the whole interface correspond 
to changes in the amount of vapor near the hydropho- 
bic plate, / ~ .001 passes -1 corresponds to the lifetime 
of vapor fluctuations, similar to the conclusion of Zhang 
et.al. The fact that the relaxation time of fluctuations 
of the entire interface corresponds to the frequency at 
which power-law behavior begins in S(f) suggests that 
interface fluctuations, rather than spontaneously- forming 
vapor bubbles, are the cause of the fluctuations in c(t). 

If D is the smallest length scale of fluctuations in the 
interface relevant to fluctuations in the contact density, 
then contributions to h(x,y;t) from Fourier modes with 
wavelengths smaller than D should be minimal. A new 
function, h£>{x,y\t), with no contributions from short- 
wavelength Fourier modes, can be created by taking the 
inverse transform of equation JSJ, but cutting off the sum 
over wavevectors at ko = 1/D: 

ko ko 

kx=0 ky=0 

When the small-wavelength Fourier modes are not in- 
cluded, hD{x,y;t) is no longer restricted to take on dis- 
crete values. The lack of high-frequency Fourier modes 
also smoothes out the interface sufficiently that it seldom 
comes in contact with the hydrophobic plate, so the def- 
inition of c(t) in equation @ is no longer useful in the 
context of fluctuations of hi){x,y;t). Instead, a clipping 
function is used to create a function cd(£), 

c D (t) = J2e[\-h D (x,y;t)] (8) 

x,y 

where is the Heaviside function, and the smooth in- 
terface h[)(x, y \ t) is considered to be in contact with the 
hydrophobic plate when < A. A is chosen for 

each D so that the time average of Cd (t) equals the time 
average of the total contact density c(t). 

To estimate the smallest- wavelength fluctuations of the 
interface relevant to the fluctuations in the c(t) time se- 
ries that was collected during the simulation, the func- 
tions h D (t) with D = 6,7,8,9, 10, 12, 14, and 16 a were 
computed from the time series of h(x, y; t) collected dur- 
ing the simulation, and the corresponding A's ranged 
from A = 1.50 for D = 6 a, to A = 1.67 for D = 16 a. 
The power spectra So(f) = |c.d(/)| were computed, and 
the slope of each curve was measured on a log-log plot 



between / = 0.001 passes -1 and / = 0.015 passes -1 , 
to see how many Fourier modes could be removed from 
h(x,y;t) but still obtain a slope statistically indistin- 
guishable from the slope of S(f) over the same frequency 
range. It should be noted that this calculation was car- 
ried out as a technique to analyze the data collected dur- 
ing the simulation, and not to make a statement about 
the exact power-law behavior of S(f), since the data are 
quite noisy. The value of the slope for each D is dis- 
played in Tabled The largest value of D with a slope still 
within one standard deviation of the true slope of S(f), 
— 1.48 ±0.04, is D = 8 a, which suggests that fluctuations 
with wavelengths smaller than 8 a are not important to 
the fluctuations in c(t). 

When a similar calculation was performed where the 
small- 1 k| (large- wavelength) Fourier modes of the inter- 
face were progressively removed from the sum in equa- 
tion JSJ, the slope of the resulting power spectrum was 
far outside of the standard deviation of the true slope. 
In order to estimate the largest length scale of interface 
fluctuations relevant to fluctuations in c(t) , a short simu- 
lation was run on a larger lattice. A simulation identical 
to that described in Section [n] but with L = 1024, was 
run for approximately 2000 Monte Carlo passes through 
the lattice. The power spectrum of the interface fluctu- 
ations, \h(k x , k y ; t)\ 2 was calculated every 10 passes, and 
the time-averaged power spectrum is shown in Figure 
plotted versus k — |k|. The power begins to depend on k 
at approximately k — .01 a , which suggests that fluc- 
tuations over length scales up to approximately 100 a are 
relevant in the simulation. Thus, it seems that interface 
fluctuations over length scales between 8 a (~ 1.5 nm) 
and 100 a (~ 19 nm) cause the fluctuations in c{t) in the 
simulation. 

Additional simulations were performed to determine 
the effects of the parameters of the model upon the re- 
sults described above. Simulations at higher and lower 
temperatures (corresponding to physical temperatures 
between 290 and 320 K), and with different values of 
(3a a, were carried out. (3e and (3[i were not varied in- 
dependently so that the system would remain close to 
liquid-vapor coexistence, and I3<jb was not varied because 
the interface is so far from the hydrophilic plate that the 
results would be unchanged. 

The results of these parameter changes were consistent 
with physical intuition. As temperature increases, the 
interface fluctuates more, causing larger-magnitude fluc- 
tuations in c(t) and also causing {h) to increase (that is, 
the interface moves farther from the hydrophobic plate). 
The power in the spatial Fourier modes of the interface 
increases, and the leveling off of the power spectrum at 
small k moves to smaller fc, indicating that larger length- 
scale interface fluctuations are important, as would be 
expected. The behavior of S(f) does not change notice- 
ably over the temperature range studied, and the linear 
regime begins at approximately the same value of /. 

/3 a a was varied from 0.0 to 0.2/3e. Increasing (3(ja 
corresponds to making the hydrophobic plate more hy- 
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TABLE I: Slope of the linear regime on a log-log plot of Sd (/) 
for different values of D (the units of D are lattice spacings 
a). The slopes have a standard deviation of approximately 
0.04. 



D slope 

6 -1.50 

7 -1.48 

8 -1.46 

9 -1.42 

10 -1.41 
12 -1.40 
14 -1.36 
16 -1.34 



drophilic, and as would be expected, doing so decreases 
(h) and the magnitude of the fluctuations in c(t), and in- 
creases (c) . The power in the interface Fourier modes de- 
creases substantially, and the power spectrum levels off at 
much larger values of k, indicating that large- wavelength 
Fourier modes become less important in the interface dy- 
namics, which is reasonable since the interface is so close 
to the plate and thus has little room in which to move. 
When j3(TA is increased to 0.2/3e, (h) = 0.89, indicating 
that the hydrophobic plate is almost completely wetted, 
and the linear regime in S(f) disappears. 



interface, and support the hypothesis that the large fluc- 
tuations and power spectrum of viscous response that 
were seen experimentally can be attributed to the fluc- 
tuating liquid- vapor interface interacting with competing 
hydrophobic and hydrophilic walls. The power spectrum 
of fluctuations in c(t) has apparent power-law behavior 
at intermediate frequencies, but the power law is a less 
rapidly-decaying function of frequency then the apparent 
experimental value of /~ 2 . However, the combination of 
the noise in the power spectrum from the simulations, 
and the lack of a precise experimental value, means there 
is no firm basis for a detailed discussion of the exact value 
of the exponent. Regardless of that exact value, a Fourier 
analysis of the liquid- vapor interface suggests that fluctu- 
ations of the interface over length scales between 1.5 and 
19 nm cause the fluctuations in c(i) that were observed 
in the simulation. It is probable that fluctuations over 
these length scales would be averaged out in a simulation 
of a system with the much larger plates studied experi- 
mentally, but it is interesting that such large fluctuations 
could be observed in such a simple system. 
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FIG. 1: Cross section of the lattice in the simulation. A pos- 
sible configuration of vapor (white) cells and liquid (shaded) 
cells is also shown. 
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FIG. 2: A portion of the time series of the contact density, 
normalized by the time-averaged value (c) = 2673. The value 
of the average contact density indicates that approximately 
25% of cells in the z = 1 layer have m = 1 at any given time. 
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FIG. 3: The power spectrum S(f) resulting from the time 
Fourier transform of c(t). 
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FIG. 4: The time autocorrelation function for fluctuations 
in the k = (0, 0) spatial Fourier mode of the liquid-vapor 
interface. 
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FIG. 5: Time averaged power spectrum of interface fluctua- 
tions as a function of k = |k|. 



